Geometric Heat Flux for Classical Thermal Transport in Interacting Open Systems 
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We study classical heat conduction in a dissipative open system composed of interacting oscillators. By 
exactly solving a twisted Fokker-Planck equation which describes the full counting statistics of heat flux flow- 
ing through the system, we identify the geometric-phase-like effect and examine its impact on the classical 
heat transport. Particularly, we find that the nonlinear interaction as well as the closely related temperature- 
dependence of system-parameters are crucial in manifesting the geometric-phase contribution of heat flux. Fi- 
nally, we propose an electronic experiment based on RC circuits to verify our theoretical predictions. 
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Understanding the general features of transports is one of 
the main goals in non-equilibrium statistical physics. Among 
many others, time-dependent driven transports like driven par- 
ticle (mass, probability) transport and driven heat conduction 
are attracting an increasing attention. In particular, the latter 
is of special interests [1-9] because of both its theoretical and 
practical importance in phononics [10], where one may uti- 
lize temporal modulations to alternatively achieve flexible dy- 
namic control of thermal energy in various phononic devices 
[10]. 

In driven quantum systems, an intrinsic geometric contribu- 
tion of the phase of a wave function will be emergent [11] and 
it has been proved to have profound effects on many physical 
properties including thermal related ones [12-14]. Similar ge- 
ometric contributions have been also uncovered in the driven 
transport of noninteracting particles [15-17]. These pioneer- 
ing efforts culminated in the discovery of geometric-phase- 
like contribution in generating functions [18, 19], which then 
inspired the identification of geometric heat flux in a single 
quantum junction [9]. Namely, even under slow modulations, 
heat flux and fluctuations are not merely a simple temporal 
average of their static counterparts, but contain an extra geo- 
metric contribution regardless of driving rates. 

Unlike noninteracting particle transport, heat conduction in 
solid is typically modeled by interacting lattices with reser- 
voirs, where energy is transported in the absence of particle 
flow [20]. Therefore, whether and how the geometric con- 
tribution can emerge in classical heat conduction is still an 
open question. Moreover, the nonlinearity (anharmonicity) 
has been found of special importance in phononic devices 
[10]. However, the role of nonlinear interaction in the man- 
ifestation of geometric heat flux is still not yet explored, al- 
though works about time-dependent classical heat conduction 
in interacting lattices have already been carried out [2-7]. 

In this Letter, we shall address the above mentioned ob- 
jectives by exactly solving a twisted Fokker-Planck equation, 
which describes the full counting statistics of heat flux flow- 
ing through a classical open system of interacting oscillators. 




FIG. 1: (color online). A sketch of two coupled classical Brownian 
oscillators in contact with two Langevin heat baths. 



We identify the geometric -phase effect on generating func- 
tions and examine its impact on the classical heat transport. In 
particular, we find that the nonlinearity of interaction as well 
as the related temperature-dependence of system parameters 
are crucial to the manifestation of geometric heat flux. Other- 
wise, for a linear system without temperature-dependent pa- 
rameters, the geometric -phase effect is absent or only observ- 
able in high order heat fluctuations. Furthermore, by pointing 
out the analogy of a coupled RC electric circuit and interact- 
ing oscillators, we are able to implement an electric experi- 
ment to verify the theoretical predictions of geometric-phase 
effects in heat transport. 

We start with a typical interacting open system: two cou- 
pled Brownian oscillators in contact with two heat baths, as 
shown in Fig. 1. The vibrational dynamics is described by a 
set of Langevin equations: rriiXi = —d Xi V(xi,X2) — Jiii + 
£i, (i = 1, 2), where mi and Xi are the mass and displacement 
of oscillator i. V(xi,X2) denotes the interaction potential. 7; 
depicts the viscosity of bath i, or say, the coupling strength 
between oscillator i and bath i. is the white noise with vari- 
ance (Cdt^jit 1 )) = 2'y i T i 6 ij 6(t - t'), where T t is the tem- 
perature of bath i. Q(t) = J d xl V{x\, x-i)i\dt' is defined to 
describe the heat transferred from heat bath 1 to 2 during time 
t. This vibrational dynamics is a typical interacting model for 
heat transport which has been also used to describe the well- 
known Feynman ratchet-pawl model [21]. For microdynam- 
ics at low Reynolds number where oscillators either possess 
nano-sizes or move in an extremely viscous media, the inertia 
terms are negligible [22], and in turn we have an overdamped 
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dynamics: 



Ti^i + dx 1 V(x 1 ,x 2 ) = Ci, 
72*2 + d X2 V(xi,x 2 ) = £>• 



(1) 



For an harmonic coupling V(xi,x 2 ) = k{x\ — x 2 ) 2 /2 with 
spring constant k, the system can be described by a Fokker- 
Planck equation [23]: 

d P (y,t) [ (T x T 2 \ d 2 ( k k\ d 1 
dt Y\li 72 / \7i 12/ oy J 

(2) 

where = .Ti(i) — i2(i). Then the transferred heat from 
bath 1 to 2 within time t is Q(t) = f Q dt'ky(^i — ky)/-yi. 

To study the time-dependent driven heat transport, we in- 
troduce the characteristic function of the joint probability 
p(y,Q,t), defined as z(y,x,t) = f^dQe 1 *® p(y,Q,t). 
This characteristic function satisfies a twisted Fokker-Planck 
equation: 

dtz(y,x,t) = L x (t)z(y,x,t), (3) 
where, following derivations from [24, 25], one can find, 



L = ( Tl I T2 ) d 2 , I k ( 1 I 1 2Tl i\\v d 
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The twisted Fokker-Planck operator L x (t) is generally time- 
dependent, wherein parameters k(t), 7^ (i) and Tj (t) could be 
subject to adiabatically cyclic modulations. The "adiabatic" 
here means the modulation period T p is much larger than the 
system's characteristic time scale of relaxation, T c , namely 
[26]: 



T p » T c = 

My 
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(5) 



It is easy to verify that Eqs. (3) and (2) have the same 
initial condition z(y,x, 0) = p(y,0). When \ — 0, we 
have z(y,0,t) — p(y,t) from the definition and in turn 
Eq. (3) reduces to Eq. (2). Integrating over the degree of 
freedom y in z(y,x,t), we obtain the characteristic func- 
tion of Q: Z( X ,t) = Jdyz(y,x,t) = J dQe ixQ P(Q,t), 
with P(Q,t) = J dyp(y,Q,t). Thus, the cumulant 
generating function is G(\) = linif^oo t^ 1 In Z(x, t), 
which generates the n-order cumulant of heat fluctuations: 

Iim t _>oo«O n ))A = a? x G(x)lx=o- 

Following [18], the cumulant generating function of adi- 
abatically driven system can be separated into two parts - 
the dynamic contribution and the geometric one: G(\) = 
limt-^oc \ \i\Z{x, t) = G d yn + Ggcom- The dynamic contri- 
bution, Gdyn = sp- Jq" dt\ (\, t), with A (x, t) denoting the 
ground-state eigenvalue of L x (t), survives whenever system 
parameters are static, or experience single or multiple modu- 
lations; whereas the appearance of the geometric contribution 
Ggcom requires at least two parameter modulations. For the 



case of periodically driving pairs (ui(t), u 2 (t)), which could 
be chosen from any two of k, -fj , Tj, we have [26] 



Ggcom = —=r jl duxdu^u^ix), 
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where ^o(V'o) denotes the corresponding left(right) ground- 
state eigenfunction, and the subscript uiu 2 denotes the in- 
tegral area enclosed by the modulating contour. Clearly, 
J r „ ltl2 (x) has the physical meaning of curvature of the pa- 
rameter space (ui, u 2 ) for L x . It is of pure geometric origin 
and is independent of the modulation speed. Mathematically, 
Ggeom is an analog of the geometric Berry phase [19] in quan- 
tum mechanics, where the wave function will gain an extra 
phase after a cyclic evolution [27, 28]. Similarly, in the full 
counting statistics of cyclic driven systems, the cumulant gen- 
erating function (analog of phase) in the exponent of the char- 
acteristic function (analog of wave function) will also gain an 
additional term [18, 19]. Both extra terms share the similar 
geometric origin from the nontrivial curvature in the system's 
parameter space. In turn the nth cumulant of heat fluctuations 
has two separate contributions as well [9, 18]: 
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Of prime interest is the first cumulant, the average heat flux 

J = Jdyn + Jgcom, with 

1 f Tp 
Jdyn = — I dt 

J seom " 
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Apparently, to study the geometric heat flux as well as 
other interest transport properties, we need to solve the eigen- 
problem of L x , which however spans the infinite-dimensional 
Hilbert space and is usually difficult to tackle. Fortunately, af- 
ter some algebra, we can cast Eqs. (3, 4) to the Schrodinger's 
eigen-problem of quantum harmonic oscillator [26], with the 
help of which we analytically obtain the ground-state eigen- 
value for L x (t), 



Ao(x,«)=r(l-0)/2, 
and the right ground-state eigenfunction, 



(11) 



r + r6 — 2i\kD x . 
Vo{y,X^) = exp ( — y ), (12) 



as well as the corresponding left eigenfunction, 



vo(y>x,t) 
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where 9 = y/1 - 4k?D 1 D 2 /r 3 (~x + «(1/T a - 1/T X )) X , 
r = fx + r 2 , D = £>i + L> 2 with r_, = fc/7j and = Tj/jj. 

Substituting Eq. (11) into Eq. (9), we see that the dynamic 
heat flux is just the temporal average of its static counter- 
part: J dyn = T~ l J Q P dtJ st (t), with J st = d ix \ \ x=0 = 
k(Ti — 2a)/(7i + 72), where the integral concerns all possi- 
ble time-dependent parameters. It is interesting to notice that 
generally given T\(t) = T 2 (t), not only the average, but all 
the odd order cumulants of the dynamic flux will vanish as 
well, due to the even symmetry Xq(x) = ^o(—X + *(l/^2 — 
1/TO) - Ao(-x). 

Substituting Eqs. (12, 13) into Eq. (7) and in turn Eq. (10), 
we can then study the geometric heat flux J g0O m under any 
pair-parameter manipulations. We first consider a special case 
of two system-bath couplings 71(f), 72(f) being modulated. 
In this case, we find the curvature = 0, that is, no mat- 

ter how arbitrarily one drives these two couplings, the geomet- 
ric contributions to all cumulants of heat transport are always 
zero. The geometric effect is absent in the case of merely 
modulating system-bath couplings. Similar phenomenon is 
also observed in quantum heat transport [9], which implies 
that this may have connection with some universal pumping 
restrictions [29] in open systems. Furthermore, this absence 
of geometric heat flux may relate to a no-pumping theorem 
of a different quantity, probability current, in closed driven 
systems without explicitly connecting with heat baths [30]. It 
shows that, in terms of our Eq. (2), the probability current is 
absent when (Tij-yi+T^J-yl) * s tinxe-independent. In fact, when 
T\ = T-2, this ratio is indeed independent of 7, so that there 
is no-pumping of probability current no matter how we drive 
(71,72)- 

In the case of modulating any other combination of two sys- 
tem parameters, however, the geometric contribution emerges. 
For example, if we modulate one system-bath coupling 72 (f ) 
and one bath temperature Ti(f), the nonzero first derivative 
of curvature, d ix T^ Tl | x =o = ~7x72/[2(7i + 72) 3 ] ± can 
induce nonzero geometric heat flux J gC om 7^ 0. For another 
example, if two baths are kept isothermal T\ = T2 = To, 
then the dynamic heat flux will be always zero, Jd yn = 0. 
Even so, we still can realize the nonzero heat transport through 
the geometric contribution by modulating, like (k, 72), so that 
J = -/gcom = 77 1 JJ km dkd l2 { 7 lT /[2fc(7i + 72 ) 2 ]}. 

The most typical modulation is driving bath temperatures 
(T 2 ,Xi) [3-6], because implementing this protocol is much 
easier than modulating k(t) and 7j(f) in practice. The curva- 
ture for this two-temperature driving reads: 



T 2 T! 



7172(72 ~ 7i)fa) 2 
2(71 + 72) 3 # 3 



(14) 



When the couplings are symmetric, 72 = 71, one can see 
that J-t 2 t x = 0. That is, the geometric -phase effect is ab- 
sent in the symmetric linear (harmonic) system. Interestingly, 
even when two couplings are asymmetric, 72 ^ 71, we find 
dix-^T-iTx |x=o — in spite of nonzero Tt 2 t x , so that the geo- 
metric heat flux is still vanishing, J gC om = 9i x G gconl \ x= o = 



0. In other words, although the geometric effect exists in the 
classical asymmetric linear system, its contribution is not ob- 
servable if one measures only the average flux. The geometric 
effect can only manifest itself in the higher order heat fluctu- 
ations, like the shot noise of currents, i.e., d? Gg eom \ x= o oc 

Recalling the modulation (72, Xi) can induce the nonzero 
geometric heat flux, we speculate that as long as the viscos- 
ity 72 is temperature dependent, the modulation (T2, Ti) can 
be effectively equal to the modulation (72, I\). As a conse- 
quence, in such systems of temperature-dependent viscosity, 
merely modulating (T2,7i) may produce nonzero curvature 
thus nonzero geometric heat flux. In fact, viscosity is gen- 
erally temperature-dependent [31]. Phenomenologically, we 
can assume jj = 70 + aT/\ with j = 1, 2. Thus for the typ- 
ical modulation (T 2 ,T\), we indeed obtain the nonzero geo- 
metric heat flux J, 



gcom 
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(15) 



Clearly, the temperature-dependence plays a key role for the 
manifestation of geometric effect. Given a modulation cycle, 
the geometric heat flux J g0 om vanishes when the temperature- 
dependency a — > 0; while J gC0 m increases to a saturated value 
as a increases. When n = 0, jj also becomes temperature 
independent, and in turn J goom becomes zero as well. 

Considering the intrinsic effective temperature- 
dependencies of system parameters are ubiquitous in 
nonlinear interacting oscillators [32, 33], we thus speculate 
that with the help of nonlinearity, the existing geomet- 
ric effect is able to manifest itself into the geometric 
heat flux as well. Therefore, we further consider the 
FPU-/3 model [10], with nonlinear interacting potential: 
V(xi,x 2 ) = ^-{xi - x 2 ) 2 + ^{xi - x 2 ) A . As we shall see 
in follows, that the nonlinear strength k 2 can induce effective 
temperature dependence of system parameters [32, 33] is the 
crucial ingredient to manifest geometric heat flux. 

In Fig. 2, we numerically simulate Eq. (1) and calculate the 
average geometric heat Q p = J Eeom T p , defined for a driving 
cycle. Given the temperature driving protocol: T 2 = 0.09 + 
0.06cos(27rf/Tp+7r/4) ) T 1 = 0.09+0.06 sin(27rf/T p +7r/4), 
we have a zero dynamic flux zero but a nonzero Q p . When T p 
becomes large (adiabatic limit), Q p saturates to a fixed value 
independent of T p , which indicates that it is purely a geomet- 
ric property. The geometric heat per driving cycle does not 
rely on the driving rate, but only depends on the geometry of 
the driving contour in parameter spaces. The deviations of Q p 
from the fixed value at the fast driving regime are due to the 
breakdown of the adiabatic precondition. 

The inset of Fig. 2 verifies that when the nonlinear inter- 
action reduces to harmonic coupling (k 2 = 0), the geomet- 
ric contribution disappears (J gcom = 0). Only with the help 
of nonlinearity, the geometric effect of temperature modula- 
tions can manifest itself into the heat flux. Moreover, in- 
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FIG. 2: (color online). Nonlinear effect on manifestation of geomet- 
ric effect in classical heat transport for FPU-/3 model, fei = 0.5, 
7i = 72 = 5. The physical value correspondences of those dimen- 
sionless units can be found in Ref. [10]. The error bar denotes the 
standard derivation of 25 times simulations, each is averaged over 
10 e periods. 

creasing k 2 can enhance Q p , until to a saturated value, which 
coincides qualitatively with the behavior of the analytic re- 
sult Eq. (15), by increasing a. In fact, from the viewpoint of 
nonequilibrium Green's functions [34], the nonlinear interac- 
tion effect in thermal transport is reflected in the temperature- 
dependent effective self-energies, which in our case are ex- 
actly the temperature-dependent ji . 

Although we focus on a two-coupled-oscillator system at 
the moment, it could be straightforward to generalize the 
above analysis into arbitrarily long coupled-oscillator model 
with inertial terms, of which the eigenvalues and eigenvec- 
tors of the twisted Fokker-Planck operator can be obtained in 
terms of appropriate phonon Green's functions [25]. 

Now, we shall examine the previous studies to see whether 
it is justified for neglecting geometric heat flux. In [3, 4], 
only one bath temperature is under cyclic manipulation so 
that the absence of geometric heat flux is justified, because 
for cyclic driving, the manifestation of geometric effect re- 
quires at least two parameter modulations in order to enclose 
a nonzero area in the parameter space. The same is true 
in [7], where the interacting lattice is only cyclically driven 
by one mechanical force. In [4—6], although two bath tem- 
peratures are modulated, they are varying either isothermally 
Tx = T 2 = T + AT(t) [5] or inversely T x = T a + AT(t), 
T 2 = T - AT(t) [4, 6]. In this way, the modulation con- 
tour only closes as a line segment. Therefore, the absence of 
geometric heat flux is also justified. 

In the last part, we would like to propose an experimental 
implementation to demonstrate our predictions on manifes- 
tation of geometric effect in classical heat transport through 
coupled oscillators. In view of the well-known electric anal- 
ogy of interacting oscillators' Brownian motion [23, 35], we 
are able to map the oscillator system into a RC circuit, where 
two resistors of resistance Rj are arranged in parallel with a 
capacitor of capacitance C, as shown in Fig. 3(a). The left 
(right) circuit part is subjected to a thermal reservoirs of tem- 
perature Ti (T 2 ), which generates a Gaussian voltage fluctu- 
ation SVi (5V2), so called Johnson-Nyquist noise, with vari- 
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FIG. 3: (color online), (a) A parallel RC electric circuit subject to 
Johnson-Nyquist noises. The analogy with two coupled Brownian 
oscillators with heat baths is reflected by the parameter correspon- 
dences: (qj, Rj, 1/C, SVj) o (xj, 7j, k, £.,). (b) Manifestation 
of geometric effect for classical heat transport in RC electric cir- 
cuit with temperature-dependent resistances, (c) Q p — J gC0 mTp as 
a function of AT. The simulations are well fitted by ttAT 2 /(8T ), 
approximated from Eq. (17). Parameters are a = lOkfi/K, C = lpF, 
To = 300K. 

ance (SVi(t)SVj(t')) = 2R i T i 5 ij S(t~t') [23, 35]. q 3 denotes 
the charge going through the resistor Rj and dqj/dt is the 
corresponding electric current. Consequently, the dynamics 
of this RC circuit is described by: 

R 1 dq 1 /dt+(q 1 -q 2 )/C = 6V 1 , 

(16) 

R 2 dq 2 /dt+(q 2 -q 1 )/C = 6V 2 . 

Note that Eq. (16) is of a similar form as the overdamped 
dynamics Eq. (1) for interacting oscillators. The heat trans- 
ferred over time t now is analogously defined as Q(t) = 
Jo <7i (?i ~~ 12)/Cdt' , which actually is the charging work done 
by the left reservoir on the capacitor. Thereafter, we can mod- 
ulate C(t), Rj(t),Tj(t) to test our predictions. 

Nonlinear capacitor can be used to simulate the nonlinearity 
effect. To mimic the related temperature-dependent viscosity, 
we can choose resistors of temperature-dependent resistance, 
e.g. Rj = aTj, with contacting them to their respective reser- 
voirs. From Eq. (15), we then have 

Iff T T 

Jgcom"// dT 2 dT x — — — 3 ■ (17) 

For cyclic modulation: T 2 = T a + AT cos{2nt/T p + tt/4), 
T = To + ATsin(27rt/T p + tt/4), we will have zero J dyn 
but nonzero J g0 om- Although this geometric heat flux is in- 
dependent of C and Rj, the adiabatic precondition Eq. (5) 
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requires T p 3> T c = CRiR^/(Ri + #2). Assume we use 
a = 10 kfi/K, C = 1 pF, T = 300 K, we will have the 
adiabatic condition T p 3> T c « 3 /is. Analytic results are 
verified by numerical simulations of Eq. (16), as plotted in 
Fig. 3(b). Clearly, for T p = 20 /is, the geometric heat per 
cycle already reaches the adiabatic limit. Therefore, as long 
as the system evolves a long time, i.e. with large number of 
cycles, we can accumulate a large geometric heat, e.g., for 
T p = 20 fj,s, AT = 50 K, Q p = 0.283 meV, after each minute 
we can have 847 eV. Moreover, if we scale up many RC cir- 
cuits in parallel, we are able to obtain even larger geometric 
heat flux. Alternatively, by increasing the modulation ampli- 
tude, one can also increase geometric heat flux, as shown in 
Fig. 3(c). 

Given the fact that the fluctuation theorems in electric cir- 
cuits [36] as well as the energy rectification have been exper- 
imentally realized in a nonlinear electrical transmission line 
[37], we believe that our prediction of geometric -phase effect 
on heat transport can be also experimentally validated in a 
foreseeable future. Furthermore, as our ability to design and 
manipulate nano/micro-sized systems improves, we believe 
that the present study of geometric energy (heat) flux could 
provide a new means of energy harvesting by harnessing the 
ubiquitous cyclic changes in the universe. 
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SUPPLEMENTARY MATERIAL FOR "GEOMETRIC HEAT FLUX OF CLASSICAL THERMAL TRANSPORT IN 

INTERACTING OPEN SYSTEMS" 



In this supplement, we are going to (1) analytically solve the eigen-problem of the twisted Fokker-Planck equation, (2) expose 
the condition of so-called "adiabatic", (3) detail the derivations of geometric-phase effect in generating functions, which finally 
leads to the geometric heat flux. 
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Exact solution of the twisted Fokker-Planck equation 

Transport behaviors in the long time limit are of our central interest. They are governed by the ground state of the twisted 
Fokker-Planck operator L x . We thus need first exactly solve the eigen-problem of the twisted Fokker-Planck equation, with 
time-independent L x : 



d t z = L x z, 



with 



L x = k B 



(-- 






V7i 


72 / 





1 1 2k B T 1 . 



7i 72 



7i 



d 



1 



Note that in the main text, we set k B = 1 for the sake of clarity. We make an ansatz of the solution: 

°° v- 



n=0 



1 



k B Ti , . ,2 

tx)y +k 

7i 7i / V7i 72 



(18) 



1 k B T x 



7i 



-IX 



(19) 



(20) 



where C n is the coefficient depending on initial conditions, A n denotes the n-th eigenvalue of L x , e~ «" f n (y) is the correspond- 
ing eigenfunction, and T is a parameter to be determined in follows. Put this ansatz into the twisted Fokker-Planck equation, 
and utilize the following relations: 



d_ 

dy 
d 2 

dy 



-fn(y) 



i^2 e ~ 4T My) 



2T MV) + e Qy , 



1 



y 



-gye w f n {y) + -^e ^ /„(„) - — e « 



y „=£df n (y) y _-£.df n (y) 



dy 



2T 



dy 



dy 2 



(21) 
(22) 



we then have 



k B 



7i 72/ dy 2 
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7i 
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To eliminate the term related to yd y , we found that one needs to set 

T 

Thus, the above equation reduces to: 

k 2 ( i + i 

\7l 72 





fa 




\7l 




J_ 




72 



2k B T 1 
7i 



(24) 



k B 
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V7i 
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— + 

dy 2 



2k B Ti 
71 



4fc s (£ + f ) 

Let us set r± 2 = fc/71 2> -D1.2 = k B T 12 /ji^ then we have: 

(r!+r 2 -2ixfc^i) 2 



2 / fcsTi 2 1 . 

fc I (zx) IX 

7i 7i 



If— + — ) 

2 \7i 72/ 

(25) 



+ (ixf k 2 D 1 - zxfcri 



4(D 1 + L» 2 ) 



y2 + n+12 _ K fn[y) = 0j 



(26) 



_c^_ (n + r 2 ) 2 
dy 2 4 (£>! + D a ) 



Akix 



(ri + r 2 )' 



(D 2 n - Dir 2 - ixkD x B 2 ) 



n + r 2 



K}My) = 0- (27) 



Then set r = r x + r 2 , D = D 1 +D 2 ,6= Jl + ^^^(i(l/T 2 - 1/T 2 ) - x), we can simplify the above equation as 



9 2 



(28) 



We further set Y = ^y and F n (Y) = f n (y), which finally leads to the eigen-problem of quantum harmonic oscillator' 
Schrodinger's equation: 



(29) 



This equation requires | (l — ^r 2 -) = 1 + 2n, with n = 0, 1, 2, From any textbook of quantum mechanics which solves the 

Schrodinger's equation of quantum harmonic oscillator, we know the eigenvalue: 



and the eigenfunction 



A„ = -[!-(! + 2n) 



where H n denotes the n-th order Hermite polynomials. 

Since we already have T = D/(r — 2i\kDi), as a consequence, the ansatz of the solution now reads: 



n=0 



Substitute it back to the twisted Fokker- Planck equation 3 t z = L x z, we have 



rt 

Yd' 



r(l+9)-a» x fcgl -.2 . / rt) . , r(l+e)-2» x fcgl ..2 

ivC * D v H n U—y) = \ n e " D 



[r0 
V V 2L> ' 



such that A n is indeed the eigenvalue of the operator L x , the right eigenfunction for the operator L y 



is. 



^n(y,X)=e ™ y H n {^j—y). 

Straightforwardly, the corresponding left eigenfunction, which is bi-orthonormal with the right one, reads, 

1 



tpn(y,x) 



rQ r(i-e)-2i X kp 1 2 rr I r6 , 



2 n n\ V 2ttD 

This left eigenfunction of L x , is just the right eigenfunction of the adjoint operator L+ which has the property [1]: 

L^(p„(y,x) = K<p n (y,x), 

corresponding to 

L x ip n (y,x) = Kipn(y,x), 

so that we have the scalar product 



f + oo 



f+oa 



dy<p m {y,x)L x ip n (y,x) = dyL+(p m (y,x)tpn{y,x) = dy^ n (y 7 x)L x f m {y,x)- 



We also may thus normalize the functions according to 

r+oo 



dy<p m {y,x)^n{y,x) = 5 n 



(30) 



(31) 



(32) 



(33) 



(34) 



(35) 



(36) 



(37) 



(38) 



(39) 
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Under the adjoint operation "+" [5], one can easily find that 



d d d d . d 

ay ay ay ay ay 



Therefore, from Eq. (19), we have 



T+ , T x T 2 \ d 2 ,/l 1 2k B Tx. \ d , 2 fk B Ti ,. , a 1. \ 2 ,k B Ti. 
\7i 72/ ay z \7i 72 7i / cty V 7i 7i / 7i 

Then, following the similar procedure, one is able to arrive at the left eigenfunction Eq. (35). 

Behaviors in the long time limit are of our central interest. They are governed by the ground state of the twisted Fokker-Planck 
operator L x , of which the eigenvalue, Ao(x), possesses the least negative real part. In other words, for time-independent L x , 
lim t _ i . 00 Z(x,t) ~ e'* 10 ^'* and in turn lirnt_ >00 ((Q n ))/i = 9" x Xq(x)\x=o- If tne system is under time-dependent modulation 
such that the Fokker-Planck operator is time-dependent L x (t), then as discussed in the main text, the full counting statistics rely 
on the instantaneous right and left ground-state ipo(y, x, t) an d VaiUi X> t)- Therefore, we only needs the information about the 
ground state. From Eqs. (30, 34, 35), we finally have the instantaneous ground state: 




r + r6 -2ixkD 1 2 
V>o{y, X, t) = exp I — y 

Similar solutions in a Fokker-Planck equation without the counting parameter x were discussed in [1,2]. Similar results about 
the eigenvalue and the right eigenfunction, but with the left eigenfunction absent, were given by [3]. 

Condition for "adiabatic" 

From Eq. (3), we write down the first three term for the sake of clarity: 

z(y, X ,t) = C a e x ° t + C 1 e Xlt + C 2 e x * t + ... 

= C e A «* [ 1 + £i e -( A «- Al > 4 + ^ e -(Ao-A 2 ) t + J 5 (42) 
\ Co Co J 

where C n are some unimportant coefficients here. Clearly, l/(Ao(x) — ^n(x)) depicts the characteristic relaxation time of 
the n-th mode to the ground state. Let us set the counting parameter x = to look into spectrum of the operator L x= q in 
the real physical space. From Eq. (30), we can see that the spectrum is ordered by n and the eigenvalue of the ground state 
Ao(x = 0) = 0. This is apparent, since it corresponds to the long time steady state of the evolution. So, the rate of relaxation 
to the steady state is determined by the inverse energy gap l/(Ao(x = 0) ~ M(x = 0)). As a consequence, the system's 
characteristic time scale of relaxation is given as 

7172 (43) 



Ao(x) ~ Ai(x) 



k(ji + 72) ' 



where Xi(x — 0) = fc/71 + fc/72 has the physics meaning of the effective damping rate. Therefore, given the system, which 
is already in its steady state after a long time evolution, as long as we drive it slowly such that the driving period T p 3> T c , the 
system can always dwell in its steady state (ground state). In other words, T p 3> T c is right the so-called "adiabatic" condition. 

Geometric phase contribution in cumulant generating functions 

We assume the system is already at its steady state, or say ground state. And then, time-dependent modulations are imposed 
on the system adiabatically. Thus, we can make the ansatz z(y, x, t) = Co(t)e^> Mx,* )dt ^ an( j substitute it into the 

time-dependent twisted Fokker-Planck equation dtz(y, x, t) = L x (t)z(y, x, t). We then obtain 

Co (t)eti x " ( ^ dt ' Mv,X,t) + C (t)e$ x ° (x ' 4 ' )dt ' i> a (y, X , *) = 0. (44) 
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Utilizing the bi-orthonormal condition dyip m (y 7 Xi t)ipn(y> X: t) 
that: 



S mn > we left multiply <po, and do integral over y, such 



C (t) = -C (t) / dy(p (y,x,t)ipo(y,x,t)- 

J —oo 

Therefore, in the adiabatic limit, we obtain 

Co(t) = C (0)e-^ dt 'fZo d yMy,x,t')4o(y,x,t')_ 

Consequently, the characteristic function finally reads 



z(x,t) 



dyC o (0)Mv,X,t) 



exp 


i 







X (x,t')dt' 



exp 



- / dt' 



dy<po(y,x,^)iPo(y,x>^) 



(45) 



(46) 



(47) 



The first exponent, the time integral of the instantaneous ground-state eigenvalue, is an analog of the dynamic phase. While the 
second additional exponent resulting from the time-evolving of ground eignstates is an analog of the geometric phase. Different 
from the conventional phase of the wave function in quantum mechanics, here the "phase" refers to the cumulant generating 
function in the exponent of the characteristic function, which will contribute to the full counting statistics of the quantities of 
interest. Similar geometric phase contribution to the generating function so as to the full counting statistics is firstly discovered 
by Sinitsyn and Nemenman in discrete chemical kinetics [4]. 

Successively, the cumulant generating function of adiabatically driven systems can be separated into two parts - the dynamic 
phase contribution and the geometric phase contribution: G(x) = lirnt-^oo j In Z(x, t) = Gd yn + G gcom , with 



G 



dyn 



Ga 



1 

— / dtx ( x ,t), 

1 p Jo 
1 



T, 



dt 



p Jo 



dypo(y,x,t)ipo(y,x,t). 



(48) 
(49) 



Note that here the contribution of initial conditions tr 1 In [ dyCo{0)ipo(y, x, t)] is assumed negligible in the long time limit. 
The dynamic phase contribution survives whenever system parameters are static, or experience single or multiple modulations. 
While the existence of the geometric phase contribution requires at least two parameter modulations. For the case of periodically 
driving (ui(t),U2(t)), which could be chosen from k, jj, Tj, using Stokes theorem, we have 



G& 



1 



duidu 2 j r UlU2 (x), 



where the subscript U1U2 denotes the integral area enclosed by the modulating contour of (iti(i), ^(i)), and 

'dip dip dipodtpo' 



,(X) 



dy 



du\ du2 dv,2 du\ 



(50) 



(51) 



is a classical analog of the quantum mechanical Berry curvature. Different from the curvature usually defined for discrete Hilbert 
space of matrix-like Hamiltonian, the curvature here is in the continuous function space, indicated by the integral over infinity. 
Clearly, the curvature F UlU2 (x) is also of pure geometric origin, since it is independent of the modulation speed. 
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